Does interfertility decrease with phylogenetic distance in Schiedea?

Binary pollen fertility data

Data from Weller et al. 2001. F1 hybrids that made greater or less than 50 percent fertile pollen.

fert <- read.table(file="data/fertility.csv", sep="\t", header=T)
fert$sp1 <- toupper(fert$sp1)
fert$sp2 <- toupper(fert$sp2)
fert$g50 <- fert$g50+1
fspecies <- data.frame(name=levels(factor(c(fert$sp1,fert$sp2))))
g <- graph_from_data_frame(fert, directed=F, vertices=fspecies)
g_edges = get.edgelist(g)

Clustering

Phylogeny on the right, clustering by network of fertile pollen on left

#clustering by network
adj <- get.adjacency(g, attr="g50")
adj[adj==2] <- 0.5
g.d <- graph_from_adjacency_matrix(adj, weighted=T, mode="undirected")
ceb <- cluster_edge_betweenness(g.d)

#trim phlyo
sch.f <- sch
sch.f$tip.label <- toupper(paste("",substr(sch.f$tip.label,1,4),sep=""))
sch.f<-drop.tip(sch.f,sch.f$tip.label[-match(fspecies$name, sch$tip.label)])

#cophenetic plot
library(phytools)
ceb_tree <- cophylo(as_phylo(ceb), sch.f)
Rotating nodes to optimize matching...
Done.
plot(ceb_tree)

Arcplot

F1 hybrids that made greater (green) or less than (red) 50 percent fertile pollen.

#phylo-arcplot
par(mfrow=c(1,2))
par(mar=c(1.2,0,1.2,0))
plot(sch.f)
par(mar=c(0,0,0,0))
arcplot(g_edges, show.labels=F, horizontal=F, ordering=sch.f$tip.label, col.arcs=fert$g50+1)

Fertility vs. phylogenetic distance

adj <- as.matrix(get.adjacency(g, attr="g50"))
adj[upper.tri(adj, diag=T)]<-NA
adj[adj==0] <- NA
phydist <- cophenetic.phylo(sch.f)
phydist[upper.tri(phydist, diag=T)]<-NA
pf <- na.omit(data.frame(pd=c(phydist), ft=c(adj-1)))
cor(pf$ft, pf$pd, use="complete.obs")
[1] -0.3740546
cor.test(pf$ft, pf$pd, use="complete.obs")

    Pearson's product-moment correlation

data:  pf$ft and pf$pd
t = -2.9085, df = 52, p-value = 0.005331
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.5833851 -0.1181252
sample estimates:
       cor 
-0.3740546 

Logistic regression

#Logistic regression 
model <- glm(ft~pd, family=binomial(link='logit'), data=pf)
summary(model)

Call:
glm(formula = ft ~ pd, family = binomial(link = "logit"), data = pf)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)    4.546      1.689   2.691  0.00712 **
pd          -177.869     70.976  -2.506  0.01221 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 71.188  on 53  degrees of freedom
Residual deviance: 62.457  on 52  degrees of freedom
AIC: 66.457

Number of Fisher Scoring iterations: 4
confint(model)
                  2.5 %     97.5 %
(Intercept)    1.660361   8.395113
pd          -337.364722 -54.308760
wald.test(b = coef(model), Sigma = vcov(model), Terms = 2)
Wald test:
----------

Chi-squared test:
X2 = 6.3, df = 1, P(> X2) = 0.012

Logistic regression plot

plot <- ggplot(pf, aes(pd, ft)) + stat_smooth(method="glm", method.args=list(family="binomial"), se=T, color="black") + scale_y_continuous(breaks=c(0, 1), labels=c("<50%", ">50%")) + scale_x_continuous() + xlab("Phylogenetic distance") + ylab("Hybrid pollen fertility") +  scale_fill_gradientn(colors=c("white","black")) + theme_bw()

plot + geom_bin2d()

Continuous pollen fertility data

crosses <- read.table("data/hpol_crosses.csv", sep="\t", header=T, colClasses = list(cross="factor",crosstype="factor"))
str(crosses)
'data.frame':   120 obs. of  11 variables:
 $ cross    : Factor w/ 120 levels "108","11","116",..: 13 24 105 107 118 2 4 9 20 23 ...
 $ mompop   : int  847 847 847 847 847 847 847 847 847 847 ...
 $ momid    : chr  "B" "B" "B" "B" ...
 $ dadpop   : int  794 794 794 794 794 842 842 842 847 791 ...
 $ dadid    : chr  "3" "3" "3" "3" ...
 $ momsp    : chr  "adam" "adam" "adam" "adam" ...
 $ dadsp    : chr  "hook" "hook" "hook" "hook" ...
 $ mombs    : chr  "d" "d" "d" "d" ...
 $ dadbs    : chr  "h" "h" "h" "h" ...
 $ crosstype: Factor w/ 4 levels "1","2","3","4": 3 3 3 3 3 4 4 4 4 4 ...
 $ notes    : chr  "" "" "" "" ...
pol <- read.table("data/hpol.csv", sep="\t", header=T, colClasses = list(cross="factor", id="factor",crosstype="factor"))
pol$vp <- pol$viable/200
pol$mombs <- factor(ifelse(pol$crosstype %in% c("1","2"), "H", ifelse(pol$crosstype=="z", "z", "D")))
pol$dadbs <- factor(ifelse(pol$crosstype %in% c("1","3"), "H", ifelse(pol$crosstype=="z", "z", "D")))
pol$momgrp <- with(pol, factor(paste(momclade,momisle,mombs)))
pol$dadgrp <- with(pol, factor(paste(dadclade,dadisle,dadbs)))
pol$momsp <- crosses$momsp[match(pol$cross, crosses$cross)]
pol$dadsp <- crosses$dadsp[match(pol$cross, crosses$cross)]
str(pol)
'data.frame':   2239 obs. of  16 variables:
 $ cross    : Factor w/ 113 levels "108","11","116",..: 99 99 99 99 99 10 10 10 10 10 ...
 $ id       : Factor w/ 38 levels "1","10","11",..: 25 25 25 25 25 11 11 11 11 11 ...
 $ viable   : int  122 165 151 166 173 189 191 183 184 179 ...
 $ crosstype: Factor w/ 5 levels "-9","1","2","3",..: 5 5 5 5 5 5 5 5 5 5 ...
 $ date     : int  1 1 1 1 1 -9 -9 -9 -9 -9 ...
 $ momisle  : chr  "O" "O" "O" "O" ...
 $ dadisle  : chr  "O" "O" "O" "O" ...
 $ momclade : chr  "2" "2" "2" "2" ...
 $ dadclade : chr  "4" "4" "4" "4" ...
 $ vp       : num  0.61 0.825 0.755 0.83 0.865 0.945 0.955 0.915 0.92 0.895 ...
 $ mombs    : Factor w/ 2 levels "D","H": 1 1 1 1 1 1 1 1 1 1 ...
 $ dadbs    : Factor w/ 2 levels "D","H": 1 1 1 1 1 1 1 1 1 1 ...
 $ momgrp   : Factor w/ 13 levels "1 K H","1 N H",..: 5 5 5 5 5 11 11 11 11 11 ...
 $ dadgrp   : Factor w/ 16 levels "1 K H","1 N H",..: 13 13 13 13 13 5 5 5 5 5 ...
 $ momsp    : chr  "adam" "adam" "adam" "adam" ...
 $ dadsp    : chr  "keal" "keal" "keal" "keal" ...
heatmap(with(pol, table(momsp,dadsp)), col=inferno(512), scale="none", Colv=NA, Rowv=NA)

#polmeans <- aggregate(vp~cross*momgrp*momclade*momisle*mombs*dadgrp*dadclade*dadisle*dadbs, data=pol, FUN=function(x) c(mean(x), sd(x), length(x)))
polmeans <- aggregate(vp~momsp*dadsp, data=pol, FUN=function(x) c(median(x), quantile(x,0.25), quantile(x,0.75), length(x))) #mean is median!!!
polmeans <- do.call(data.frame, polmeans)
colnames(polmeans)[3:6] <- c("mean", "q25", "q75", "N")
#write.csv(polmeans,  "results/hpol_means.csv", row.names=F)
#alllvl <- sort(unique(c(levels(polmeans$momgrp), levels(polmeans$dadgrp))))
#grps <- cbind(alllvl, read.table(text=alllvl, sep=" ", col.names=c("clade","isle","bs")))
species <- unique(c(pol$momsp,pol$dadsp))
grps <- data.frame(sp=species) 

library(reshape2)
poltbl <- dcast(polmeans, momsp ~ dadsp, value.var="mean", fun.aggregate = mean)
rownames(poltbl) <- poltbl$momsp
poltbl$momsp <- NULL
poltbl2 <- as.matrix(poltbl)
poltbl2[is.nan(poltbl2)] <- 0
ggplot(pol, aes(x=momsp, color=dadsp, y=vp)) + geom_boxplot()

ggplot(pol, aes(x=momclade, color=dadclade, y=vp)) + geom_boxplot() + geom_point()

ggplot(pol, aes(x=momisle, color=dadisle, y=vp)) + geom_boxplot() + geom_point()

ggplot(pol, aes(x=momclade, color=dadclade, y=vp)) + geom_boxplot() + geom_point()

ggplot(pol, aes(x=reorder(cross, viable, FUN=mean), y=vp, color=crosstype)) + geom_boxplot() + geom_point()

ggplot(pol, aes(x=crosstype, y=vp, color=dadclade, fill=momclade)) + geom_boxplot()

#heatmap(poltbl2, col=viridis(512), scale="none", Rowv = NA, Colv=NA)

#ggplot(polmeans, aes(x=factor(momsp, levels=alllvl), y=factor(dadgrp, levels=alllvl), fill=mean)) + geom_tile() + scale_x_discrete(drop=F) +  scale_y_discrete(drop=F) + scale_fill_viridis_c()
#ggplot(polmeans, aes(x=momclade, y=dadclade, fill=mean)) + geom_tile() + scale_fill_viridis_c()
ggplot(polmeans, aes(x=momsp, y=dadsp, fill=mean)) + geom_tile() + scale_fill_viridis_c()

hist(pol$viable, breaks=40)

colfunc <- viridisLite::plasma

g <- graph_from_data_frame(polmeans[,c("momsp","dadsp","mean","q25","q75","N")], directed=T, vertices=grps)
lgw <-  layout_in_circle(g)
#plot.igraph(g, layout=layout_in_circle(g, order=order(grps$clade)), edge.color=colfunc(100)[round(E(g)$mean*100)],  vertex.color=rainbow(10)[grps$clade], vertex.label.color=rainbow(10)[as.integer(grps$isle)+6], vertex.label.cex=2)
plot.igraph(g, layout=layout_in_circle(g))

gu <- as.undirected(g, mode="collapse", edge.attr.comb = "mean")
adj <- get.adjacency(gu, attr="mean")
g.w <- graph_from_adjacency_matrix(adj, weighted=T, mode="undirected")
#plot(g.w, layout_with_fr(g.w), edge.color=colfunc(100)[round(E(g)$mean*100)], edge.width=2,   vertex.color=rainbow(10)[grps$clade], vertex.label.color=rainbow(10)[as.integer(grps$isle)+6], vertex.label.cex=2)
plot(g.w, layout=layout_with_fr(g.w), edge.color=colfunc(100)[round(E(g)$mean*100)], edge.width=2)

#dendPlot(cle, mode="phylo", use.edge.length = T, main="leading eigen")
plot(as.dendrogram(cluster_leading_eigen(gu)))

g_edges = get.edgelist(g)
arcplot(g_edges, show.labels=T, horizontal=F, col.arcs=colfunc(100)[round(polmeans$mean*100)], sorted=T)

Clustering

ceb <- cluster_leading_eigen(gu)
#trim phlyo
sch.f <- sch
sch.f$tip.label <- tolower(paste("",substr(sch.f$tip.label,1,4),sep=""))
sch.f<-drop.tip(sch.f,sch.f$tip.label[-match(species, sch.f$tip.label)])
species1f <- species1
rownames(species1f) <- tolower(rownames(species1f))
species1f <- species1f[sch.f$tip.label,]
#BScols <- setNames(c("black", "navy", "lightblue", "white", "white", "white"),1:6)
BScols <- setNames(RColorBrewer::brewer.pal(6,"Set2")[c(6,4,2,5,1,3)],1:6)
#intraspecies <- setNames(polmeans$mean[match(paste(sch.f$tip.label,sch.f$tip.label), paste(polmeans$momsp,polmeans$dadsp))], sch.f$tip.label)
intraspecies <- setNames(polmeans$mean[match(paste(names(V(g)),names(V(g))), paste(polmeans$momsp,polmeans$dadsp))], names(V(g)))


#cophenetic plot
ceb_tree <- cophylo(as_phylo(ceb), sch.f)
Rotating nodes to optimize matching...
Done.
plot(ceb_tree)

#phylo-arcplot
#png("plots/arcplot.png", height=1500, width=3000, pointsize=56)
par(mfrow=c(1,3) , bg="grey30")
par(mar=c(1.5,0,1.5,0))
plot(sch.f, align.tip.label=T, tip.color=BScols[as.integer(species1f$Breeding.System)], edge.width=8, edge.color="grey80", cex=1.5, lwd=8)
add.scale.bar(y=0.002, x=0.002, lcol="white",col="grey80")
par(mar=c(0,0,0,0))
arcplot(g_edges, show.labels=F, horizontal=F, ordering=sch.f$tip.label, col.arcs=plasma(512, direction=-1)[round(polmeans$mean*512)], sorted=T, show.nodes=T, lwd.arcs=6, col.nodes = plasma(512, direction=-1)[round(intraspecies*512)], pch.nodes = 19, cex.nodes=1)
plot.new()
#legend(x=-0.05,y=1, legend = rev(levels(species1f$Breeding.System)), col = rev(BScols), pch = 19, cex=1.2, pt.cex=2.5, bty="n", xpd=T, text.col= rev(BScols))
library(plotrix)
color.legend(0,0.55,1,0.65,legend=seq(0,1,0.2),rect.col=plasma(512, direction=-1),cex=1, col="grey80", align="rt")
text(0.2, 0.7, "Pollen viability", col="grey80", cex=1.5)
#dev.off()

Fertility vs. phylogenetic distance

phydist <- cophenetic.phylo(sch.f)
#phydist[upper.tri(phydist, diag=T)]<-NA

adj <- as.matrix(as_adj(g, attr="mean"))
#diag(adj) <- NA
#adj[upper.tri(adj, diag=T)]<-NA
adj[adj==0] <- NA
adj <- adj[sch.f$tip.label,sch.f$tip.label]

adj.N <- as.matrix(as_adj(g, attr="N"))
adj.N <- adj.N[sch.f$tip.label,sch.f$tip.label]

adj.q25 <- as.matrix(as_adj(g, attr="q25"))
adj.q25 <- adj.q25[sch.f$tip.label,sch.f$tip.label]
adj.q75 <- as.matrix(as_adj(g, attr="q75"))
adj.q75 <- adj.q75[sch.f$tip.label,sch.f$tip.label]

comb <-function(x,y) {mapply(function(q,r) paste(sort(c(q,r)), collapse="-"), x, y)}
labels <- outer(sch.f$tip.label, sch.f$tip.label, comb)

pf <- na.omit(data.frame(pd=c(phydist), ft=c(adj), N=c(adj.N), q25=c(adj.q25), q75=c(adj.q75),crossname= c(outer(sch.f$tip.label, sch.f$tip.label, paste, sep=" x ")), crossname.sort=c(labels)))
pf$q75[pf$q75==1]<- NA # safe for logit transform
cor(pf$ft, pf$pd, use="complete.obs")
[1] -0.7844253
cor.test(pf$ft, pf$pd, use="complete.obs")

    Pearson's product-moment correlation

data:  pf$ft and pf$pd
t = -11.026, df = 76, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.8573058 -0.6807202
sample estimates:
       cor 
-0.7844253 
pdft <- ggplot(pf, aes(x=pd, y=ft, weight=200*N)) + 
    geom_line(aes(group=crossname.sort), col="grey70", linewidth=2, alpha=0.5) + 
  theme_grey() + 
  xlab("Phylogenetic distance (mutations/site)") + 
  ylab("Pollen viability")

pdft + geom_point() + geom_linerange(aes(ymin=q25, ymax=q75)) + ylim(c(0,1))

pdft + scale_y_continuous(trans="logit", breaks=c(0.01,seq(0,1,by=0.05), 0.995)) + geom_text(aes(label=crossname), size=3) # trans_breaks(function(p) { log(p/(1-p)) }, function(x) { 1/(1+exp(-x)) },  n=10)

pdft  + geom_smooth(method="glm", method.args=list(family="binomial"), color="grey70", fill="grey80") + geom_text(aes(label=crossname, color=N)) + coord_trans(y="logit") + scale_y_continuous(breaks=c(0.01,seq(0,1,by=0.05), 0.994)) + scale_color_gradient(trans="sqrt", low="darkblue", high="green2")

#scale_y_continuous(breaks=seq(0,1, by=0.05), limits=c(0,1), expand=c(0,0))

#ggsave("pdft.png", height=9, width=9)

Arcplot

#phylo-arcplot
#TODO: arrange arcs front-to-back by phydist?
library(gridBase)
library(grid)

#png("plots/arcplot.png", height=1500, width=4000, pointsize=20, res=300)
par(mfrow=c(1,4) , bg="grey30")
par(mar=c(1,0,1,0))
plot(sch.f, align.tip.label=T, tip.color=BScols[as.integer(species1f$Breeding.System)], edge.width=3, edge.color="grey80", cex=1, lwd=4, label.offset=0.0001)
add.scale.bar(y=0.6, x=0, lcol="grey80",col="grey80", cex=1, lwd=3)
par(mar=c(0,0,0,0))
arcplot(g_edges, show.labels=F, horizontal=F, ordering=sch.f$tip.label, col.arcs=plasma(512, direction=-1)[round(polmeans$mean*512)], sorted=T, show.nodes=T, lwd.arcs=3, col.nodes = plasma(512, direction=-1)[round(intraspecies[order(names(intraspecies))]*512)], pch.nodes = 19, cex.nodes=1)
plot.new()
#legend(x=-0.05,y=1, legend = rev(levels(species1f$Breeding.System)), col = rev(BScols), pch = NA, cex=1, pt.cex=2.5, bty="n", xpd=T, text.col= rev(BScols))
library(plotrix)
color.legend(0.1,0.05,0.2,0.55,legend=seq(0,1,0.2),rect.col=plasma(511, direction=-1),cex=0.45, col="grey80", align="rt", gradient="y")
text(0.26, 0.6, "Pollen viability", col="grey80", cex=1)
plot.new()    
vps <- baseViewports()
pushViewport(vps$figure) 
vp1 <-plotViewport(c(0,0,0,0)) 
print(pdft  + 
        geom_smooth(method="glm", method.args=list(family="binomial"), color="grey70", fill="grey80") +
        geom_linerange(aes(ymin=q25, ymax=q75)) +
        geom_point(aes(color=ft), size=2) + 
        coord_trans(y="logit") + 
        scale_y_continuous(breaks=c(0.01,seq(0,1,by=0.1), 0.99)) + 
        scale_color_viridis_c(option="plasma", direction=-1) + 
        theme_classic() + 
        theme(plot.background = element_rect(color = "grey30", fill = "grey30"), panel.background = element_rect(fill = "grey30", color  =  NA),  legend.background = element_rect(color = NA, fill = "grey30"), axis.text = element_text(color = "grey80"), axis.title = element_text(color = "grey80")) + 
        guides(color=F), vp = vp1)

#dev.off()

Binomial regression on means

# model <- glm(ft~pd, family="quasibinomial", weights=N*200, data=pf)
# summary(model)
#plot(model)

library(glmmTMB)
betabin.tmb <- glmmTMB(round(cbind(ft*N*200,(1-ft)*N*200))~pd, data=pf, family="betabinomial")
summary(betabin.tmb)
 Family: betabinomial  ( logit )
Formula:          round(cbind(ft * N * 200, (1 - ft) * N * 200)) ~ pd
Data: pf

     AIC      BIC   logLik deviance df.resid 
  1220.9   1227.9   -607.4   1214.9       75 


Dispersion parameter for betabinomial family (): 3.63 

Conditional model:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)    2.9620     0.3353   8.833  < 2e-16 ***
pd          -134.2268    16.5402  -8.115 4.85e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Binomial regression on raw data

phydist.m <- melt(phydist)
polmeans$pd <- phydist.m$value[match(paste(polmeans$momsp,polmeans$dadsp), paste(phydist.m$Var1,phydist.m$Var2))]
pol$pd <- phydist.m$value[match(paste(pol$momsp,pol$dadsp), paste(phydist.m$Var1,phydist.m$Var2))]
pol$tot <- 200
pol$inviable <- pol$tot-pol$viable

ggplot(pol, aes(y=vp,x=pd)) + geom_point()

# model <- glm(cbind(viable, inviable)~pd, family="quasibinomial", data=pol)
# summary(model)
# #plot(model)
# library(jtools)
# effect_plot(model, pred = pd, interval = TRUE, data=pol) + geom_point(data=pol, aes(pd,vp), color="blue")
# coef(model)

library(glmmTMB)
# plot.glmmTMB <- function(x, ...) {  qplot(x$fitted, x$residuals, xlab="Fitted values", ylab="Residuals", ...) + geom_smooth()}
# 
# betabin.tmb <- glmmTMB(cbind(viable,inviable)~pd, data=pol, family="betabinomial")
# summary(betabin.tmb)
#plot(betabin.tmb)
# plot(fitted(betabin.tmb)~pol$pd)
# points(polmeans$mean~polmeans$pd, col="red")
#coef(betabin.tmb)

betabin.tmb.momdadsp <- glmmTMB(cbind(viable,inviable)~pd + (1|momsp) + (1|dadsp), data=pol, family="betabinomial")
#load("betabin.tmb.momdadsp.Rdata")
summary(betabin.tmb.momdadsp)
 Family: betabinomial  ( logit )
Formula:          cbind(viable, inviable) ~ pd + (1 | momsp) + (1 | dadsp)
Data: pol

     AIC      BIC   logLik deviance df.resid 
 20209.5  20238.1 -10099.8  20199.5     2234 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 momsp  (Intercept) 0.8521   0.9231  
 dadsp  (Intercept) 0.5991   0.7740  
Number of obs: 2239, groups:  momsp, 15; dadsp, 16

Dispersion parameter for betabinomial family (): 3.82 

Conditional model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   1.6064     0.3249   4.943 7.67e-07 ***
pd          -81.4044     4.9036 -16.601  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#plot(betabin.tmb.momdadsp)
#save(betabin.tmb.momdadsp, file="betabin.tmb.momdadsp.Rdata")

######Plots############
#Looks like the random effects are overfitting :( GLM does the best
#pol$fitted <- fitted(betabin.tmb)
pol$fitted <- fitted(betabin.tmb.momdadsp)
#pol$fitted <- fitted(model)
polmeans$fitmean <- aggregate(fitted~momsp*dadsp, data=pol, FUN=mean)$fitted
  ggplot(polmeans, aes(pd, fitmean)) +    
    geom_smooth(aes(y=mean), method="glm", method.args=list(family="binomial"), color="red", fill="pink") +
    geom_smooth(method="glm", method.args=list(family="binomial"), color="black") + 
    geom_point() + geom_point(aes(y=mean), color="red") + 
    geom_segment(aes(xend=pd, yend=mean))  + geom_text(aes(y=mean, label=paste(momsp,dadsp,sep=" x ")), color="blue")

####Random effects#####
momsp.ef <- data.frame(sp=rownames(ranef(betabin.tmb.momdadsp)$cond$momsp), ef=ranef(betabin.tmb.momdadsp)$cond$momsp[[1]])#, sd=ranef(betabin.tmb.momdadsp, sd=T)$cond$momsp[[1]])
momsp.ef$sp <- reorder(momsp.ef$sp, momsp.ef$ef)
ggplot(momsp.ef, aes(sp,ef)) + geom_point() +# geom_linerange(aes(ymin=ef-sd, ymax=ef+sd)) + 
  ggtitle("Maternal random effects") + xlab("Species") + ylab("Effect")

dadsp.ef <- data.frame(sp=rownames(ranef(betabin.tmb.momdadsp)$cond$dadsp), ef=ranef(betabin.tmb.momdadsp)$cond$dadsp[[1]])#, sd=ranef(betabin.tmb.momdadsp, sd=T)$cond$dadsp[[1]])
dadsp.ef$sp <- reorder(dadsp.ef$sp, dadsp.ef$ef)
ggplot(dadsp.ef, aes(sp,ef)) + geom_point() +# geom_linerange(aes(ymin=ef-sd, ymax=ef+sd)) + 
  ggtitle("Paternal random effects") + xlab("Species") + ylab("Effect")